#This file takes two csv data files to create figures 1 through 3
#The first csv file is statewide covid data by day
#The second file contains estimates of experiences and anxiety generated in the stata .do file "wide analysis". See that file for generation of estimates. 

library(foreign)

data<-read.csv("Replication data for figure 1.csv")
jpeg("Figure1.jpg", width = 750, height=750, pointsize=18, quality = 100)
plot(data$day, data$sevendaymean, type = 'l', ylim=c(0,1), xlim=c(-7,449), lwd=2, cex=1.25,
	xlab="", ylab = "New Cases Per Thousand", main = "", axes = FALSE)
axis(side = 1, at=c(7,38,68,99,129,160,191,221,252,282,313,343,372,403,433),  
	las = 2, cex.axis = 1, labels = c("Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec","Jan",
	"Feb","Mar","Apr","May")) 
axis(side =2, at = seq(0, 1, .1), las = 1, 
	cex.axis = 1, labels =seq(0,1,.1))
rect(130,0,141,1, col="grey")
rect(177,0,193,1, col="grey")
rect(243,0,253,1, col="grey")
rect(303,0,309,1, col="grey")
rect(408,0,428,1, col="grey")
lines(data$day, data$sevendaymean, lwd=2)
dev.off()

data<-read.csv("Replication data for figures 2 and 3.csv")
jpeg("Figure2.jpg", width = 750, height=1500, pointsize=18, quality = 100)
par(mfrow=c(2,1))
plot(data$day, data$knowpos, ylim=c(0,.6), xlim=c(-7,449), pch=19, cex=1.25,
	xlab="", ylab = "Proportion", main = "A. Health Threats", axes = FALSE)
axis(side = 1, at=c(7,38,68,99,129,160,191,221,252,282,313,343,372,403,433),  
	las = 2, cex.axis = 1, labels = c("Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec","Jan",
	"Feb","Mar","Apr","May")) 
axis(side =2, at = seq(0, .6, .1), las = 1, 
	cex.axis = 1, labels =seq(0,.6,.1))
box()
points(data$day, data$testpos, pch=17, cex=1.25)
points(data$day, data$knowdeath, pch=15, cex=1.25)
lines(data$day, data$testpos, lwd=2)
lines(data$day, data$knowdeath, lwd=2)
lines(data$day, data$knowpos, lwd=2)
segments(data$day,data$kplb,data$day,data$kpub,lty=3, lwd=2)
segments(data$day,data$kdlb,data$day,data$kdub,lty=3, lwd=2)
segments(data$day,data$tplb,data$day,data$tpub,lty=3, lwd=2)
legend("topright", legend = c("Know someone who tested positive", 
	"Know someone who died", "Tested positive"), pch=c(19,15,17), cex=1)
plot(data$day-2.5, data$incloss, ylim=c(0,.6), xlim=c(-7,449), pch=19, cex=1.25,
	xlab="", ylab = "Proportion", main = "B. Economic Threats", axes = FALSE)
axis(side = 1, at=c(7,38,68,99,129,160,191,221,252,282,313,343,372,403,433),  
	las = 2, cex.axis = 1, labels = c("Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec","Jan",
	"Feb","Mar","Apr","May")) 
axis(side =2, at = seq(0, .6, .1), las = 1, 
	cex.axis = 1, labels =seq(0,.6,.1))
box()
points(data$day+2.5, data$hhjobloss, pch=17, cex=1.25)
lines(data$day, data$incloss, lwd=2)
lines(data$day, data$hhjobloss, lwd=2)
segments(data$day-2.5,data$illb,data$day-2.5,data$ilub,lty=3, lwd=2)
segments(data$day+2.5,data$jllb,data$day+2.5,data$jlub,lty=3, lwd=2)
legend("topright", legend = c("Decrease in income", 
	"Somoene in HH lost job"), pch=c(19,17), cex=1)
dev.off()

jpeg("Figure3.jpg", width = 750, height=750, pointsize=18, quality = 100)
par(mfrow=c(1,1))
plot(data$day-2.5, data$hltanx, ylim=c(2.5,3.5), xlim=c(-7,449), pch=19, cex=1.25,
	xlab="", ylab = "Mean", main = "", axes = FALSE)
axis(side = 1, at=c(7,38,68,99,129,160,191,221,252,282,313,343,372,403,433),  
	las = 2, cex.axis = 1, labels = c("Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec","Jan",
	"Feb","Mar","Apr","May")) 
axis(side =2, at = seq(2.5, 3.5, .1), las = 1, 
	cex.axis = 1, labels =seq(2.5,3.5,.1))
box()
points(data$day+2.5, data$ecanx, pch=17, cex=1.25)
lines(data$day-2.5, data$hltanx, lwd=2)
lines(data$day+2.5, data$ecanx, lwd=2)
segments(data$day-2.5,data$halb,data$day-2.5,data$haub,lty=3, lwd=2)
segments(data$day+2.5,data$ealb,data$day+2.5,data$eaub,lty=3, lwd=2)
legend("topright", legend = c("Health anxiety", 
	"Economic anxiety"), pch=c(19,17), cex=1)
dev.off()